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PACS 87. 16.D — Membranes, bilayers, and vesicles 

PACS 83 . 80 . Lz - Rheology: Physiological materials (e.g. blood, collagen, etc.) 

PACS 87. 17.Aa - Theory and modeling; computer simulation 

Abstract - The dynamics and rheology of suspensions of fluid vesicles or red blood cells is 
investigated by a combination of molecular dynamics and mesoscale hydrodynamics simulations 
in two dimensions. The vesicle suspension is confined between two no-slip walls, which are driven 
externally to generate a shear flow with shear rate 7. The flow behavior is studied as a function 
of 7, the volume fraction of vesicles, and the viscosity contrast between inside and outside fluids. 
Results are obtained for the encounter and interactions of two vesicles, the intrinsic viscosity of 
the suspension, and the cell-free layer near the walls. 



Introduction. — Suspensions of mesoscale particles 
in viscous liquids arc ubiquitous, with examples in biolog- 
ical systems (blood flow), home products (paints), food 
products (emulsions), and industrial processing (pastes). 
The suspended particles can be spheres, rods, fibers, flex- 
ible and semiflexible macromolecules, droplets, capsules, 
vesicles and cells. While the dynamics of rigid particles in 
suspension and their rheological behavior have been inves- 
tigated in considerable detail and are by now reasonably 
well understood [1] , much less is known about the dynam- 
ics and rheology of deformable particles, in particular in 
the semi-dilute regime, where hydrodynamic and steric in- 
teractions between the particles become important. 

The dynamics of soft objects, in particular under flow, 
depends on the physical origin of their deformability, like 
the surface tension at constant volume for droplets, the 
membrane bending rigidity at fixed volume and surface 
area for vesicles, and in addition the membrane shear elas- 
ticity for capsules and cells. Therefore, these systems have 
to be investigated independently to understand the rela- 
tion between the elasticity of the particles and the rheo- 
logical behavior of their suspensions. 

In the dilute regime, the vesicle dynamics shows tank- 
treading (TT), tumbling (TU) and vacillating-breathing 



dynamics, depending on shear rate 7 and viscosity con- 
trast A [IHS]. For TT quasi-spherical vesicles in three di- 
mensions (3D), the viscosity of a dilute suspension has 
been predicted to be pifTU] 
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as a function of excess area A = 47r[^(|pr)2/3 _ 1] aj^^i 
viscosity contrast A ~ 'qm/'Hout, where A and V are the 
surface and volume of the vesicle, r]in and rjout are the 
fluid viscosities of the inner and outer fluids, respectively, 
and (f) is the vesicle volume fraction. Thus, the intrin- 
sic viscosity 77/ = (77 — rjout) / (rjoutip) is predicted to be a 
decreasing function of A and A. Furthermore, rji is fore- 
seen to have a cusp-like minimum at the tank-treading to 
tumbling (or tank-treading to vacillating-breathing) tran- 
sition, and then to increase again with increasing A pillO] . 
This latter behavior has been also found in the numerical 
calculations of a two-dimensional vesicle by the boundary- 
integral approach [llj . 

These theoretical predictions have been tested experi- 
mentally [T2 l [T3 ] . While a decrease of 777 with increasing A 
was found in ref. |12| , in good agreement with the theoret- 
ical prediction ([T]) , in contrast an increase of 77/ was found 
in ref. [13]. However, the available experimental results 
are not conclusive for several reasons. First, vesicle sizes 
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Fig. 1: Configurations at consecutive times 7i = 424, 440, 
456, 472 (from top to bottom) of vesicles with viscosity contrast 
A = 1.0, reduced area A* — 0.8, reduced shear rate 7* = 2.0, 
and concentration = 0.14. One vesicle is colored blue for 
better visualisation of its evolution during tank-treading. See 
also movie SI for A = 2.0, A* = 0.8, 7* = 2.0, and (j) = 0.28. 



in suspensions are typically polydisperse. Second, viscos- 
ity measurements require a minimum volume fraction of 
vesicles, typically 5% to 10%, and are therefore difficult to 
extrapolate to the dilute limit [13]. Indeed, experiments 
have been performed recently [M] which demonstrate that 
vesicle interactions become relevant for the viscosity for 
around 10%. 

Therefore, we study here the rhcology of vesicle sus- 
pensions in the "semi-dilute" regime, vsfhere particle inter- 
actions are important, but particles are not yet densely 
packed into a glassy state. Our results are obtained from 
mesoscale hydrodynamics simulations of two-dimensional 
(2D) model systems, which allow the study of larger sys- 
tem sizes and longer time scales. Compared to the the- 
ory of ref. [nmn] and the simulations of rcf. |TT], our 
model includes thermal fluctuations and has the capabil- 
ity of studying systems over a wide range of vesicle con- 
centrations. Our main results concern the dependence of 
the intrinsic viscosity on viscosity contrast, shear-thinning 
behavior, displacements and angular oscillations in two- 
vesicle collisions, and the dependence of cell-free-layer 
thickness on shear rate 7. 

Method and Model. — Each vesicle in two dimen- 
sions is modeled as a chain of Np beads of mass nip, 
connected successively in a closed ring [T3], see fig. [T] 
Neighboring beads are connected to each other by an har- 
monic potential with spring constant k^ and average bond 



length tq] this keeps the perimeter length of the mem- 
brane constant, both locally and globally. Shapes and 
fluctuations are then controlled by a bending potential 
Vfc — (k/^o) J2ii^ ~ cos/3i), where /3i is the angle between 
the two bond vectors at bead i, and k is the bending rigid- 
ity. Finally, in order to keep the vesicle area A close to the 
target value Aq, a potential Va = kA^A — Ao)^/{2rQ) is 
employed, where fc^ is the compression modulus. Differ- 
ent vesicles repel each other at short distances via a shifted 
Lennard- Jones potential, which is truncated at its mini- 
mum Tcut- Newton's equations of motions for the beads 
are integrated by using the velocity- Verlet algorithm with 
time step Atp [T6j . 

The fluid is described by multi-particle collision (MFC) 
dynamics, a particlc-bascd mesoscale simulation technique 
[T7Hl9| . The two-dimensional fluid consists of Ng point 
particles of mass to, whose positions ri{t) and velocities 
Vi(i), i = 1,2, ...,Ns, are continuous variables. The evo- 
lution occurs in discrete time intervals Atg, and proceeds 
in two consecutive steps: streaming and collision. In the 
streaming step, particles move ballistically. In the colli- 
sion step, the particles are sorted into the cells of a regu- 
lar square lattice of mesh size a; all particles within each 
cell collide and exchange momentum. We employ here a 
variant of MFC, denoted as MFC-AT-t-a, which conserves 
both linear and angular momentum locally [JOlEI] and 
keeps the temperature constant [20] . The viscosity of the 
MFC-AT-|~a fluid in two dimensions is given by 
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with I = Atsy^ksT/m the mean- free path, fc^T the ther- 
mal energy, and n the average number of particles per cell 
[22] . The system of size L^ x Ly is placed between two 
horizontal walls which slide along the x direction with ve- 
locities Uujaii and —v^aii, respectively. Feriodic boundary 
conditions are used along the x direction. A bounce-back 
rule with virtual particles ensures no-slip boundary con- 
ditions at the walls [1TJ[5^. This generates a linear flow 
proflle Vx = jy with shear rate 7 = 2vwaii/ Ly. 

To describe the fluid-membrane interaction, membrane 
beads are modeled as hard disks, see fig. [1] The radius 
Tp of the disks is chosen large enough to ensure mutual 
overlap and a complete coverage of the membrane to pre- 
vent fluid particles from crossing the membrane. Since 
it is very important to conserve linear and angular mo- 
mentum for vesicles with viscosity contrast [211 [24], we 
employ the following scattering rule between fluid parti- 
cles and membrane disks. Scattering occurs only when a 
fluid particle j and a membrane disk i overlap and move 
towards each other, so that the conditions jr^ — r^l < rp 
and (r, — r^) • (vj — Vj) < are satisfied. A second disk 
/s = i ± 1 in the same membrane, with min |r/j — rj\, is 

selected to perform a three-body collision which conserves 
linear and angular momenta [24]. The MFC collision step 
is then performed only for those fluid particles which did 
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Fig. 2: The intrinsic viscosity rji = {rj — rjout) / {rjoutfj}) as a 
function of the viscosity contrast A for reduced shear rate 7* = 
2.0, reduced area A* = 0.8, and concentrations (j) ~ 0.05 (•), 
0.09 (A), and 0.14 (*). The dashed line is the interpolation to 
the data (•). The tank-treading-to-tumbling transition occurs 
at Xc ~ 3.7 for A* = 0.8 in the KS theory [2]. 



Fig. 3: The intrinsic viscosity tji = {rj — rjout) / {'rjoutff') ^s a 
function of the reduced shear rate 7* for reduced area A* = 0.8, 
concentration — 0.28, and viscosity contrasts A = 2.0 (•) and 
5.0 (o) with L^x Ly = (18.95 x 5.79)i?o, and A = 2.0 (□) and 
5.0 (A) with L-^x Ly = (15.79 x 6.84)i?o- When not visible, 
error bars are comparable with symbols size. 



not participate in the membrane scattering, in order to 
avoid multiple collisions with the same disk in subsequent 
time steps. The fluids in the interior and exterior of the 
vesicle may differ in their particle mass to control viscosity. 
Membrane disks interact with walls via bounce-back. 

In experiments with vesicles in shear flow, inertial effects 
are negligible since the Reynolds number Re = '^pR^/rjout, 
where p = nm/a? is the fluid mass density, is typically 
very small. We express the results in dimensionless quan- 
tities, such as the reduced area A* = Aq/ttRq (where 
i?o = Lo/{2tt) is the mean vesicle radius with membrane 
length Lq) and the reduced shear rate 7* ~ jrjoutRo/K. 
We set n = 10, lout = 0.0064a with kn = loiLt\/mout/m^n 
(in the following the subscripts out/in will refer to quan- 
tities outside/inside of the vesicle) . This implies that the 
viscosity contrast is A = run/rjout — Tnm/mout- We use 
the system size L^ = 18.95i?o, Ly = 5.79Ro, mean ra- 
dius Rq = 7.6a, and Vwaii such that Re < 0.2 for all the 
cases we considered with 0.4 < 7* < 10.0. Finally, we set 
TTi.jn such that 0.1 < A < 13.0, m.p = Sniout, Np ~ 480, 
Atp = Ats/64, rp = j'o = a/10, r^ut = a, k = 6.58fc_BTi?o, 
fc^ = 4 X IQ-'^kBT, kh = 3 X lO^fcsT, and A^ such that 
0.8 < A* < 0.95. This value of k gives rise to a similar 
amplitude of undulation modes as for lipid bilaycr mem- 
branes in 3D (where h^d — lOksT). With these choices 
for kA and kh, the area and the length of the vesicle arc 
kept constant with a deviation less than 1% of the target 
values for all simulated systems. 

Results. 

Suspension Viscosity. We flrst consider dilute and 
semi-dilute monodisperse suspensions of vesicles with flxed 
reduced shear rate 7* = 2.0. Systems with Ny = 2, 4, 
or 6 vesicles are studied, corresponding to concentrations 



(p = 0.05, 0.09, and 0.14, for different viscosity contrasts 
A. A few typical vesicle conflgurations for = 0.14 and 
A = 1.0 arc displayed in fig. [TJ The suspension viscosity 77 



from the component cTxy 



of 



is calculated numerically 

the stress tensor, so that ry = (Jxy/j HI- 

The relative viscosity {rj — rjout) /flout is a linear function 
of (j) for A* = 0.8 and various values of A as predicted by 
the Einstein relation [26]. In flg. [21 the intrinsic viscos- 
ity r]i is shown as a function of A for various concentra- 
tions with A* = 0.8. An increase of 77/ with the viscosity 
contrast is observed. We do not find an indication of a 
non-monotonic behavior — as predicted theoretically by 
eq. ^ in refs. PITU] for quasi-spherical vesicles in 3D and 
obtained numerically in 2D in ref. |llj — in the explored 
range 1.0 < A < 9.0 of viscosity contrasts, although the 
dynamic behavior changes from TT to TU at intermedi- 
ate values of A. In two dimensions, the TT-to-TU tran- 
sition is predicted to occur at Ac — 3.7 for A* = 0.8 in 
the Keller-Skalak (KS) theory [2]. However, thermal vesi- 
cle undulations, which are neglected in KS theory, pro- 
duce a continuous crossover from TT to TU for bending 
rigidities around k = GAksTRo [24], with A* = 0.7 and 
7* ^ 6. Thus, our simulation results of increasing 777(A) 
are in qualitative agreement with the experimental results 
of ref. [13] for semi-dilute systems. 

We consider next the behavior of monodisperse concen- 
trated suspensions with (j) = 0.28 (12 vesicles) with re- 
duced area A* = 0.8 for viscosity contrasts A ~ 2.0, 5.0 as 
a function of the reduced shear rate 7*. Two systems of 
size L^xLy = (18.95 x 5.79)Eo and (15.79 x 6.84)i?o are 
investigated, which have the same area but the latter be- 
ing a vesicle radius Ro wider than the former. The results 
of the intrinsic viscosity 77/ are displayed in fig. [31 The 
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Fig. 4: Relative vertical displacement of the centers of mass 
AY/Ro of two scattering vesicles with respect to the relative 
horizontal distance AX/Ro for A — 1.0, A* = 0.8, and shear 
rates 7* = 2.0 (•), 5.0 (o), and 10.0 (*). 



values of rji are not affected by the system size. For both 
values of A, a significant shear-thinning behavior is found 
over more than one decade in the reduced shear rate. The 
data also show that it is difficult to reach the low shear- 
rate plateau in simulations. This is due to the importance 
of thermal motion at low shear rates, but may also be re- 
lated to the broad TT-to-TU transition in 2D [23 where 
some tumbling events already appear in the TT regime. 
This shear-thinning is mainly due to the formation of cell- 
free layers near the walls, as expected from the Fahraeus- 
Lindqvist effect [27] . The formation of cell- free layers will 
be discussed in detail below. An analysis of the effective 
viscosity in the central part of the channel, as derived from 
the local shear rate, shows that shear-thinning of the core 
region, as observed in bulk red blood cell suspensions in 
3D, both experimentally [5S] and in simulations ^\, is not 
significant in 2D in the considered concentration range. 

Vesicle Interactions. Following the experimental work 
in refs. [TBUTl] , we study the interaction between two vesi- 
cles with A* = 0.8 in the TT regime (A = 1.0). Figure 
m displays the relative vertical displacement of the centers 
of mass Ay = ycmi — 2/cm2 during scattering with respect 
to the horizontal displacement AX = Xcmi ~ Xcm2, where 
ixcmi,ycmi) and {Xcm2,ycm2) are the positions of the cen- 
ters of mass of the two vesicles, for three different shear 
rates 7* = 2.0, 5.0, and 10.0. Movie S2 illustrates the 
vesicle interaction for 7* = 2.0. Figure [5] shows four im- 
portant effects: (i) When the vesicles are released from 
their initial positions in the upper and lower halves of 
the channel, they migrate towards the center due to the 
wall-induced lift force {AX/Rq < —2); (ii) upon collision, 
the vesicles are displaced and reach a maximum in their 
vertical separation AY corresponding to the small vesicle 
diameter {AX/Rq ~ 0); (iii) immediately after the col- 
lision, the vertical displacement is larger than before the 
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Fig. 5: Dynamics of two vesicles over two subsequent inter- 
action events. Deviations A9i of the inclination angles of both 
the two vesicles {i = 1,2, indicated by the two types of sym- 
bols) from the average stationary value are shown as a function 
of time for the run in fig. |4]with 7* = 2.0. The vertical lines de- 
note the times of the closest relative distance between the two 
vesicles. See also movie S2, which displays vesicle interactions 
in the time range 30 < -yi < 60. 



collision {AX/Rq ~ 2.5); (iv) the vesicle continue to mi- 
grate towards the center line {AX/Rq > 2.5). Different 
shear rates mainly determine the migration rate, but seem 
to have little effect on the collision process itself. Good 
agreement for the collision process is found with experi- 
mental results, which arc obtained for much wider chan- 
nels (see fig. 2 of ref. [I3|). The increased vertical sep- 
aration after scattering is in qualitative agreement with 
recent theoretical predictions for quasi-spherical vesicles 
and large inter- vesicle distances [30] . 

The interaction process can also monitored in time by 
considering the behavior of the distance d between the cen- 
ters of mass of the vesicles and the deviations A9i = Oi — Oo 
of the inclination angle 6i {i — 1,2) of the two vesicles 
from its average stationary value 0q. This latter is shown 
in fig. [5] for two consecutive scattering events for the same 
run of fig. H] with reduced shear rate 7* = 2.0. We ob- 
served a correlation of the tilt angles of the vesicles when 
the relative distance is at minimum, in agreement with the 
experimental results (compare with fig. 8 of ref. [11]). In 
fig.m the deviations A9i are shown as a function of the rel- 
ative displacement angle a, defined as the angle between 
the direction along the vesicles centers of mass and the 
flow direction, for the run in fig.|3]with 7* ~ 2.0. Data are 
averaged over four subsequent scattering processes. The 
maximum and minimum of A6i occur approximately at 
a ~ 37r/4 and a ~ 7r/4, corresponding to the compression 
and stretching directions of the shear flow field, respec- 
tively. This is again in agreement with the recent exper- 
imental results (see fig. 6 of ref. [14]), except for a small 
displacement of the minimum position. 
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Fig. 6: Deviations AS, of the inchnation angles of two inter- 
acting vesicles (i = 1, 2, indicated by the two types of symbols) 
from the average stationary value as a function of a, the angle 
between the direction connecting the vesicles centers of mass 
and the flow direction, for the run in fig. 2] with 7* = 2.0. Data 
are averaged over four subsequent interaction events. 



Cell-Free Boundary Layer. For red blood cells (RBCs) 
in capillary flow, it was first found by Fahraeus and 
Lindqvist ^7\ that cells are depleted from a layer near 
the vessel wall. It is now well understood that this 
is a consequence of the wall-induced hydrodynamic lift 
force on the cells [24]. For concentrated systems with 
(p ~ 0.28, we have measured the average thickness S of 
the cell-free boundary layers near the walls. S is defined 
as the time average of {di{t) + d2{t))/2 with (ii,2(0 = 
mini=i^...^Ar^, li{l, 2)(i), where Ny is the number of vesicles 
and /i(l,2) the closest distance of i-th vesicle membrane 
from either of the two walls. The results are reported in 
fig. [7] as a function of the reduced shear rate. The ratio 
5/Ro increases with 7* for both the values of the consid- 
ered viscosity contrasts. The data are consistent with a 
logarithmic growth of 6 with increasing 7* . With increas- 
ing system width, the values of S grow, as also observed in 
two-dimensional simulations of RBC-like vesicles in capil- 
lary flow pTj . 

The existence of boundary layers is also supported by 
considering the fluid mass density profiles in the steady 
state averaged along the flow direction x. The mass den- 
sity is lower in the cell-free boundary layers at the walls 
due to the absence of vesicles with a heavier fluid inside 
(see the inset of flg. [7]) . This effect of the boundary layer is 
also evident in steady-state velocity profiles, which display 
a smaller effective shear rate in the center and a higher 
shear rate near the wall, as compared to the imposed shear 
rate 7 (see the inset of fig. [7]). 

Finally, the behavior of (5 as a function of A for 7* = 2.0 
is shown in fig. [S] for the two system sizes. It is evident 
that there is a pronounced non-linear dependence of S on 
A, with the boundary layer decreasing at high values of 



Fig. 7: The ratio of the average thickness 5 of the vesicle-free 
boundary layers to the vesicle radius _Ro as a function of the 
reduced shear rate 7* for reduced area A* — 0.8, concentration 
(f) = 0.28, and viscosity contrasts A = 2.0 (•) and 5.0 (o) with 
L^ X Ly = (18.95 X 5.79)i?o, and A = 2.0 (□) and 5.0 (A) 
with Lx X Ly = (15.79 x 6.84)_Ro. Error bars are comparable 
with symbols size. The full line is a logarithmic fit to data 
points (•). Insets: (left) Fluid mass density and (right) fluid 
velocity component Vx, averaged along the flow (x) direction, 
both across the channel with A — 2.0 and 7* = 10.0. (Right) 
The full line corresponds to the imposed shear flow profile. 



the viscosity contrast, attaining a maximum near the TT- 
to-TU transition. This behavior can be related to the 
dependence of the lift force on the viscosity contrast, which 
has been shown [24] to be a decreasing function of A. In the 
TT phase, the increase of S might be due to the reduction 
of the tilt angle with increasing A, which facilitates the 
sliding of vesicles past each other and allows them to be 
squeezed more easily into the center of the channel. In 
the TU phase, tumbling of vesicles is suppressed near the 
wall, which further reduces the lift force |24l . 

The cell- free layer thickness S is found to grow with Ly, 
in agreement with simulations of RBCs in cylindrical mi- 
crochannels in 3D (compare flg. 10 of ref. [35]). S is found 
to increase with 7, in agreement with the results of RBC 
simulations in 3D presented in flg. 8 of ref. [35] , but at odds 
with the simulation results in fig. 11 of ref. [32] and exper- 
imental results of ref. [34] . This apparent discrepancy can 
be partially resolved by taking a closer look at the inves- 
tigated range of shear rates. In our 2D case, the increase 
of S occurs for 7* < 7; in ref. [33] , it is seen for 7* < 20 
in 3D; and in ref. [32], S is found to slightly decrease for 
7* > 4 (in refs. [32j[33] the average shear rate is used, 
computed from the average velocity in a Poiseuille flow; 
for RBCs in three dimensions, we use r = tjqRq/k with 
mean radius Ro = 3.4/xm, bending rigidity k = bOksT, 
and the plasma viscosity 770 = 0.0012 Pa s — which im- 
plies T = 0.22s — to determine the dimensionless shear 
rate). Although it is of course difficult to compare results 
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Fig. 8: The ratio of the average thickness 5 of the vesicle- 
free boundary layers to the vesicle radius Ro as a function of 
the viscosity contrast A for reduced area A* — 0.8, reduced 
shear rate 7* = 2.0, concentration = 0.28, and system sizes 
L^xLy = (18.95 X 5.79)i?o (•) and L^ x Ly = (15.79 x 6.84)i?o 
(o). Error bars are comparable with symbols size. 



of 2D and 3D systems quantitatively, we can conclude that 
there is a critical reduced shear rate 7* ~ 10 below which 
S is increasing with 7*, and above which d is constant or 
slowly decreasing. The value of 7* depends, of course, on 
the channel width and the vesicle volume fraction [32 ; the 
value above should be valid for volume fractions around 
0.3 and channel widths of about three vesicle diameters. 

Summary and Conclusions. We have investi- 

gated the dynamical behavior of semi-dilute suspensions 
of vesicles or red blood cells under shear flow in a nar- 
row gap between two walls. The advantage of the Couettc 
geometry compared to Poiseuille flow is that wall effects 
and effects of a non-linear flow profile do not interfere. We 
find the intrinsic viscosity to increase monotonically with 
increasing viscosity contrast, a pronounced shear-thinning 
behavior with increasing shear rate due to the Fahraeus- 
Lindqvist effect, displacements and angular oscillations in 
two- vesicle collisions in good agreement with experiments, 
and an increase of the cell-free-layer thickness with shear 
rate 7 below a critical reduced shear rate 7* ~ 10. 



Fruitful discussions with D. Fedosov, 1. O. Gotzc, S. 
Messlingcr, H. Noguchi, and M. Pcltomaeki are gratefully 
acknowledged. 
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